!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!! Linear solver wrapper
!!
!! \n
!!
!! This module bundles all the linear solvers, thus simplifying using
!! them and calling the respective constructors. See also the
!! comments in \c mod_linsolver_base.
!!
!! \section general General layout
!!
!! The main criteria used in this module are as follows.
!! <ul>
!!  <li> The solution of a linear system always involves three steps:
!!   <ol>
!!    <li> factor the linear system: \c factor; sometimes the
!!    factorization can be further splitted in two steps: analysis
!!    and numeric factorization, so \c factor has the following
!!    interface:
!!    <ul>
!!     <li> if the optional argument \c phase is not present, both
!!     analysis and numeric factorization are excuted
!!     <li> if the optional argument \c phase is present, then it must
!!     be either <tt>"analysis"</tt> or <tt>"factorization"</tt>;
!!    </ul>
!!    notice that after the first call to \c factor, any subsequent
!!    call wich requires an analysis must be preceded by a call to \c
!!    clean
!!    <li> solve the linear system: \c solve; this function can be
!!    called any number of times following a call to \c factor
!!    <li> clean up: \c clean; this function can be called any number
!!    of times, regardless of any previous call to \c factor and \c
!!    solve.
!!   </ol>
!!   These steps are specified as \c deffered type bound procedures of
!!   the basic type \c c_linpb, so that each solver will specify them.
!!   An additional field holds the (pointer to the) solution of the
!!   linear system, which must be a type extending \c t_stv.
!!  <li> To define a linear solver, it is necessary to privide an
!!  abstract type extending \c t_linpb. This type should implement the
!!  three subroutines \c factor, \c solve and \c clean and add all the
!!  reuired field which will allow the final user to define his/her
!!  linear problem.
!!  <li> The final user will have to define a type which extends the
!!  chosen linear solver, where he/she will provide the required
!!  information. A variable of this type will then be created to
!!  collect the unknown, the linear problem and the solver.
!!  <li> When using MPI based parallel solvers, <em>all</em> the calls
!!  to functions in the present module must be made by all the MPI
!!  processes (but the behavior depends on the chosen solver).
!! </ul>
!!
!! \section solvers Supported solvers
!!
!! Currently the following options are implemented (refer to the
!! specific modules for further details):
!! <ul>
!!  <li> <a
!!  href="http://www.cise.ufl.edu/research/sparse/umfpack/">UMFPACK</a>
!!  solvers (direct)
!!  <li> <a href="http://graal.ens-lyon.fr/MUMPS/">MUMPS</a> solvers
!!  (direct), <em>requires MPI</em>
!!  <li> the iterative solvers of <tt>mod_iterativesolvers</tt>.
!! </ul>
!!
!! \subsection ddc_solvers Parallel execution
!!
!! Some solvers support a parallel solution of a given linear system,
!! other do not. In general, when called by multiple MPI threads, the
!! behaviour of the solver routines will depend on the chosen solver.
!! Currently we have that
!! <ul>
!!  <li> UMFPACK: each thread solves locally the specified system and
!!  there is no communication.
!!  <li> MUMPS: all the threads solve in parallel a unique linear
!!  system; the system matrix cab be centralized or distributed, as
!!  specified in \c linsolver_mumps_parms and documented in the MUMPS
!!  manuals.
!!  <li> the iterative solvers of <tt>mod_iterativesolvers</tt>: all
!!  the threads solve in parallel a unique linear system.
!! </ul>
!! \warning While in serial execution all the solvers should deliver a
!! similar solution, when using multiple MPI threads a change in the
!! linear solver might imply a change in the problem solved.
!!
!! \section usage Using a solver
!!
!! The general method to define a linear problem with an associated
!! solver is to define a type which extends the chosen solver abstract
!! type and to declare a varable of the extended type. When more than
!! one linear problem must be defined, they can use different types
!! (if the solvers are different) or different variables of the same
!! type (if the linear problems are different but the solver is the
!! same).
!!
!! To be able to select the linear solver during the execution, one
!! will typically
!! <ul>
!!  <li> define a type for each solver that must be available (all
!!  these types will extend, indirectly, \c c_linpb)
!!  <li> define an <tt>allocatable</tt> variable of type \c c_linpb
!!  which will then be allocated to the type of the chosen solver.
!! </ul>
!! The main point here is that, once the variable has been allocated,
!! the calls to \c factor, \c solve and \c clean do not depend on the
!! chosen solver.
!!
!! \section details Implementation details
!!
!! The linear solver might use some pointer fields to simplify sharing
!! memory consuming data among linear problems, however the pointed
!! variables are not changed by the solvers. All the working arrays
!! are local to the solvers.
!<----------------------------------------------------------------------
module mod_linsolver

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_linsolver_base, only: &
   mod_linsolver_base_constructor, &
   mod_linsolver_base_destructor,  &
   mod_linsolver_base_initialized, &
   c_linpb

 use mod_iterativesolvers_base, only: &
   mod_iterativesolvers_base_constructor, &
   mod_iterativesolvers_base_destructor,  &
   mod_iterativesolvers_base_initialized, &
   c_itpb

 use mod_gmres, only: &
   mod_gmres_constructor, &
   mod_gmres_destructor,  &
   mod_gmres_initialized, &
   gmres

 use mod_mumpsintf, only: &
   mod_mumpsintf_constructor, &
   mod_mumpsintf_destructor,  &
   mod_mumpsintf_initialized, &
   c_mumpspb

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_linsolver_constructor, &
   mod_linsolver_destructor,  &
   mod_linsolver_initialized, &
   c_linpb, c_itpb, c_mumpspb, &
   gmres

 private

!-----------------------------------------------------------------------

! Module variables

 logical, protected ::               &
   mod_linsolver_initialized = .false.
 character(len=*), parameter :: &
   this_mod_name = 'mod_linsolver'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_linsolver_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if(mod_messages_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if( mod_linsolver_base_initialized .or. &
mod_iterativesolvers_base_initialized .or. &
                mod_gmres_initialized .or. &
            mod_mumpsintf_initialized ) then
     call error(this_sub_name,this_mod_name,                     &
       'The specific linear solve modules are already initialize')
   endif
   if(mod_linsolver_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   call        mod_linsolver_base_constructor()
   call mod_iterativesolvers_base_constructor()
   call                 mod_gmres_constructor()
   call             mod_mumpsintf_constructor()

   mod_linsolver_initialized = .true.
 end subroutine mod_linsolver_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_linsolver_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_linsolver_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   call        mod_linsolver_base_destructor()
   call mod_iterativesolvers_base_destructor()
   call                 mod_gmres_destructor()
   call             mod_mumpsintf_destructor()

   mod_linsolver_initialized = .false.
 end subroutine mod_linsolver_destructor

!-----------------------------------------------------------------------
 
end module mod_linsolver

